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ABSTRACT 

We present 3-D models of dust distribution around (3 Pictoris that produce the best fits to the 
Hubble Space Telescope Advanced Camera for Surveys' (HST/ACS) images obtained by Golimowski 
and co-workers. We allow for the presence of either one or two separate axisymmetric dust disks. 
The density models are analytical, radial two-power-laws joined smoothly at a cross-over radius with 
density exponentially decreasing away from the mid-plane of the disks. Two-disk models match the 
data best, yielding a reduced of ~1.2. Our two-disk model reproduces many of the asymmetries 
reported in the literature and suggests that it is the secondary (tilted) disk which is largely responsible 
for them. Our model suggests that the secondary disk is not constrained to the inner regions of the 
system (extending out to at least 250 AU) and that it has a slightly larger total area of dust than the 
primary, as a result of slower fall-off of density with radius and height. This surprising result raises 
many questions about the origin and dynamics of such a pair of disks. The disks overlap, but can 
coexist owing to their low optical depths and therefore long mean collision times. We find that the 
two disks have dust replenishment times on the order of 10^ yr at ~100 AU, hinting at the presence of 
planetesimals that are responsible for the production of 2"^^ generation dust. A plausible conjecture, 
which needs to be confirmed by physical modeling of the coUisional dynamics of bodies in the disks, 
is that the two observed disks are derived from underlying planetesimal disks; such disks would be 
anchored by the gravitational influence of planets located at less than 70 AU from (3 Pic that are 
themselves in slightly inclined orbits. 

Subject headings: formation: disks, planets, dynamics - dust, scattering function 
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Ow ing to its proximity (19.28 ± 0.19pc; ICrifo et ahl 
Il997f ). j3 Pictoris is one of the best studied examples 
of a Main Sequence (MS) star with a circumstcUar 
dust disk. (3 Pic was shown to have a substantial mid- 
and far-infrared excess through observ ations using the 
Infra re d Astronomical S atelite ( IRAS ) (|Aumann et al. I 
[1981 lAumann I [l985t IGillet I [HH), believed to be 
indicative of thermal radiation from ^100 K dust 
orbitin g the star. Follow-up imaging bv lSmith fc Terrili 
(|1984f ) revealed an almost edge-on system with cir- 
cumstellar nebulosity due to light scattered from the 
dust around the star. Since then, many investigators 
have mapped the surface brightne ss and color of the 
disk in optical w avelen g hts (.Smith fc Terrild [l987l; 
Artvmowicz et al.l 119891 : iLecavelier Des Etangs et al.l 



19931: iGol imowski et a lT 119931: iKalas fc JewittI 119951: 



Heap et al . 2000; Goh mowski et al.ll2006i r 



The extensive studies of the disk have also revealed 
several asymmetries in the scattered light contours pro- 
duced by t he dust surround i ng /3 Pic. These were sum- 
marized bv IKalas fc Jewit^ ()1995D and include: 1. the 
size asymmetry - the north-east (NE) extension of the 
disk stretches further out from (3 Pictoris than the south- 
west (SW) extension, 2. the surface brightness asymme- 
try - the brightness profile along the disk's spine (mid- 
plane) is broader for the NE than the SW extension, 
3. the width asymmetry - the SW extension is thicker 
than the NE extension, 4. the wing tilt asymmetry - the 
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two extensions of the disk are not perfectly aligned (i.e. 
the difference in the position angles of the midplanes of 
the opposing extensions is less than 180°) and 5. the 
butterfly asymmetry - isophote curvature is asymmetric 
across the midplane of the disk, with the asymmetry it- 
self inverted across the minor-projection axis. Several 
attempts have been made to expl ain why these asymme- 
tries arise. iKalas fc Jewittj ()1995() noted that the size and 
surface brightness asymmetries might be a direct result 
of the width asymmetry. Since the SW extension of the 
disk is more vertically extended, the SW midplane might 
appear fainter then the NE midplane and thus explain 
why the latter appears to have a larger radial size. 

Further insight into the (3 Pic system has come 
from optical and mid- infrared observations of its disk. 
iHeap et al.l (|2000) and Goli mowski et al.l (|20060 imaged 
the disk with high-resolution HST optical cameras and 
detected a warp in the inner part of the disk (^20-100 
AU from (3 Pic) indicating the presen ce of a secondary 
disk inclined at ^ 5° t o the primary (iHeap et al. 2000; 
IGolimowski et aP l2006l ) . Intriguingly, IGolimowski et aTl 
(I2006D noted that the projected spine of the secondary 
is aligned with the isophotal inflections (the butter- 
fly as ymmetry) reporte d at l arge distances from the 
star. IGolimowski et al] (|2006t ) also placed constraints 
on the brightness profiles along the spines of the pri- 
mary, secondary and composite (combined primary and 
secondary) disks. They found that the primary disk has 
a steeper midplane profile than the secondary at greater 
distances from the star and that the composite profile is 
best described with four power laws, in contrast to previ- 
ous modelling efforts that employed two power laws (eg. 
IKalas fc JewittI I1995D . To probe the regions of [3 Pic's 
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disk in accessible with optical instruments. I Wahhai et aTl 
(|20(B imaged the disk in the mid-infrared (17.9 /xm). 
They reported a set of rings embedded in f3 Pic's disk 
whose radii vary from 14 to 82 AU and whose configura- 
tion suggests the presence of multiple planets. However, 
the ri ng structure is uncertain, since IGolimowski et all 
(|2006f ) failed to observe them in scattered light. 

Strong motivation for the study of dust distribution in 
disks is spurred by the fact that such a study reveals im- 
portant clues about the presence of unseen larger bodies 
(planetesimals and planets) . For exarn ple, it was sug- 
gested that the HST warp ()Heap et al.| [2000) is caused 
by the presence of a giant planet on a slig htly inclined or- 
bit le ss than 50 AU from the primary (ex. iMouillet et al.l 
Il997f) . This prediction was possibly vindicated recently 
by the announcement of an 8Mjup candidate companion 
separated by 8 AU from (3 Pic, which, if confirmed, would 
make it one of the first planets imaged orbiting an A type 
MS star, and the planet imaged with the smallest pro - 
jected physical separation to-date (jLagrange et al.ll2^009D . 

The system's many peculiarities have yet to be ade- 
quately explained, and thus continued studies of the dust 
distribution are a necessity that might lead to the dis- 
covery of additio nal companions. Given the po ssibility 
of dust migration (iTakcuchi fc Artv mowicj|200ir ) and/or 
dust avalanches (jGrigoricva ct al. 2003), the dust and the 
planetesimal/planet distributions may not be the same. 
However, detailed modeling of micron-sized dust includ- 
ing its dynamical interaction with radiation and dust- 
dust collisions, is a prerequisite for understanding the 
distribution of larger bodies. Our work provides the first 
step for such prospective modeling, the description of the 
dust distribution. 

In this paper we presen t the best-fit two-disk m odel 
to the images obtained bv IGolimowski et al.l (|2006f) . In 
section 2 we discuss the model assumptions as well as 
the ACS images that are used for the fitting procedure. 
We dedicate section 3 to the results and section 4 to the 
discussion and analysis. 

2. MODELING AND DATA 

2.1. Multi- Parametric Model 

We use an axisymmetric multi-parameter model to de- 
scribe the distribution of the dust surrounding (3 Pic. 
Its simplicity and the fact that it has been success- 
fuUy applied in the past on more t han one occasion (eg. 
lArtvmowicz et al.lll989t iKalas fc^ewitti il995') allow one 
to quickly proceed from the numerical development stage 
to testing the model on astronomical data. The im- 
provements on the previous applications of the multi- 
parametric model include the introduction of an addi- 
tional axis of rotation (which allows us to model disks 
whose axis can point in any direction), implementation of 
a minimization algorithm (in this case a Markov chain 
Monte Carlo method) and the addition of the secondary 
disk which we fit simultaneou sly with the primary. Fur- 
thermore, recent IGolimowski et al.l (j2006[ ) observations 
have produced the most photometrically accurate maps 
to date which allow one to probe dust distribution in all 
directions (specifically at large scale heights), allowing 
one to introduce additional constraints on the model. 

In order to describe t he dust distribution, we follow 
lArtvmowicz et all ()1989f ) . We assume that both disks are 
axisymmetric. In cylindrical co-ordinates {r,z,9) then. 



the scattering cross section per unit volume is given by: 
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This form arises as a result of the vertical exponential 
drop off of the dust density. The rate of this drop off 
is determined by the parameter p. W{r) is the width 
profile: 
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This profile guarantees that the disk is thicker at 
greater radii than at smaller radii. If 7, the fiare in- 
dex, is more than 1, the disk will flare. The scale height 
Aq ensures that the disk has a specific thickness to ra- 
dius ratio at r = Rm, the location of the power law break 
describing the vertical optical thickness. 

The vertical optical thickness is proportional to the 
function T(r). r(r) has the following dependance on dis- 
tance from the centre of the star: 



T(r) = 



Vir/Rmr^"- + (r/Rm)-'^ 



(3) 



The inner radial index a and the outer radial index f3 
suppress the vertical optical thickness below and above 
the characteristic radius Rm- Only two indices are used 
since the midplane surface brightness profile of the pri- 
mary disk is known to have a single break and th us re- 
quires two power laws to describe it ([Art vmowicz et al.l 
119891: IKalas fc Jewitt|[l995HGohmowski et al.,.200d r 

To compare the brightness contours produced by our 
model with images of /3 Pic's disk, we convert the cylin- 
drical symmetric dust distribution into the observed 
brightness isophotes in the following manner. Light in- 
tensity along the line of sight through the disk is given 
by: 



I{x,y) = 



%{x,v,l)f{e)dl 



(4) 



Here x and y are the pixel co-ordinates of an image 
and 1 is the depth along the line of sight. In the above 
equation, fiux entering a unit volume having co-ordinates 
(r,z) is Lj,/[47r(r^-|-z^)], where L,y denotes the luminosity 
of P Pic at frequency v. This luminosity, along with the 
distance to the /3 Pic, is hidden inside our normalization 
parameter rj. f{9) represents the phase function and is 
described in section 2.2. The formulae for r, z and 6 
as functions of x, y and / follow from the appropriate 
co-ordinate transformations. 



[{I cosi -I- y sini)^ -I- x 
z — y cos i — I sin i 



,1/2 



cos 



-'(~l/{x^ + y^ + l'f'' 



(5) 
(6) 

(7) 



The theoretical image obtained via the described trans- 
formation will show a disk that is aligned with the major 
projection axis. Angle (j) is introduced so that this image 
can be rotated around the line-of-sight until it is aligned 
with a disk, which may not lie along the major projection 
axis (which is certainly the case for the secondary). 



Dust Distribution Around (3 Pictoris 



3 



0.1 



0.01 



S 0.001 



le-04 



: 1 1 1 


1 1 1 1 1 : 








isotropic \ 




1 1 1 


g=0.5 : 


g=0.67 ■ 
1 1 1 1 1 



20 40 60 80 100 120 140 160 180 
Phase angle (degree) 

Fig. 1. — The scattering efficiency as a function of phase an- 
gle is plotted for empirical zodiacal (Lcincrt et al. 1976) and sev- 
eral Henyey-Greenstein fits (Hcnyey & Grecnstein. ,1941i) . including 
isotropic scattering (g=0). 

2.2. Phase Functions 

In order to obtain a theoretical image, one must also 
include a phase function [f{9) in equation 4] which de- 
scribes how the light is scattered in different directions 
compared to the line-of-sight by the orbiting dust^. For 
our simulations we use d two sc attering profiles, the em- 
pirical zodi a.cal (Leinert et al.l fl976.) an d the Henyey- 
Greenstein (|Henvev fc GreensteinI Il941h . The former 
profile describes light scattered by the midsized particles 
(10-1 00 /im) with rough surfaces (jSchiffer fc ThielheimI 
|1982|) and it strongly favors forward scattering. The em- 
pirical zodiacal is proportional to: 

f{e) = 0.3(0.2 + 61/2)-^ + 1.4(61/3.3)'^ + 0.2 . (8) 

The Henyey-Greenstein scattering profile is an analyti- 
cal phase function commonly used to describe light scat- 
tered by small grains: 



1(0) 



1 



47r(l + .g2 - 25C0s6')3/2^ 



(9) 



where —1 < g < 1 is the profile asymmetry parameter. 
Positive values of g produce forward scattering, negative 
backward while <? = gives us the isotropic scattering. In 
Fig[T] we present the zodiacal light and several Henyey- 
Greenstein phase functions. In the simulations (as well 
as FiglT]) the phase functions are normalized so that their 
intergrals over the full solid angle are equal. This is done 
so that a meaningful comparison can be made between 
rj parameters obtained from different simulations. 

2.3. Data 

In order to constrain the t wo-disk model of Pic's 
dust distribution, we use the Goli mowski et al.l (|2006l ) 
ACS Hubble Space Telescope (HST) images. The ACS' 
excellent optics and HST's well studied PSF allow one 
to produce images whose quality (photometric accuracy 
and morphological details) are unparalleled by any other 
hardware combination currently in operation. Further- 
more, the ACS' wide field-of-view (FOV) [- 29" x 26"], 

^ We assume the same optical properties of dust everywhere. 



its ability to access large scale heights [~ 1.5" circu- 
lar coronographic mask only blocks out the innermost 
30 AU] and its abihty to use three filters [F435W(B), 
F606W(broad V) and F814(broad I)] to probe chromatic 
dependence introduce additional constrains on our model 
that would have been impossible with previous observa- 
tions. 

To discriminate between the features in the im- 
age associated with the disk and those intrinsic to 
the coron ographi c PSF of the disk and f3 Pic itself, 
Golimow ski et al.l (120061 ) observed the disk at two roll 
angles and obtained images of a star with similar colors 
{a Pic) prior to every exposure. Following the appro- 
priate subtractions, the combined images were decon- 
volved with synthetic PSFs produced by the Tiny Tim 
software package^ . Th roughout the reduction process 
I Golimowski et ahl ()2006f ) keep track of errors (both ran- 
dom and systematic) associated with each pixel. This 
allows them to produce meaningful error maps that in- 
clude both the systematic and statistical uncertainties in 
their data; we use these error maps for our minimiza- 
tion. The number of data points in these images once 
the coronograph (inner 30 AU) is masked out is 563136. 

2.4. Markov chain Monte Carlo 

We employ Markov chain Monte Carlo (M CMC: lFord l 
l2005f ) fitting as described for our purposes in lCroU I 20061 . 
MCMC fitting is a computationally efficient method to 
perform a full Bayesian analysis of complicted problems 
where one is required to fit a great number of often cor- 
related parameters. MCMC fitting effectively samples 
the posterior parameter distribution by performing an n- 
step intelligent random walk around the iiT-dimensional 
paramater space of interest, while recording the &t 
each point in the intelligent random walk. In addition, 
MCMC fitting allows one to explore correlations in one's 
fitting-parameters - the end result is that one is able to 
return realistic best-fit and error estimates even if the 
parameters of interest are correlated. MCM C fitting is 
thus well-suited to fit the I Golimowski et aL (2006) ob- 
servations of the (3 Pic disk with our model as it features 
a large number of parameters (i4r=18-20 parameters for 
the two-disk case) a select number of which display mod- 
est correlations. 

Flat priors in all pa rameters are used . We run our 
MCMC chains until the lGelman fc RubinI (|l99l l statistic 
is close to unity for all parameters. To determine our 
best-fit and l-tr uncertainties, w e use the ma rginalized 
likelihood method as described in lCroll I (|2006f l. 

2.4.1. Accounting for the large number of data points 

Given the large number of data-points (>500000) in 
the ACS images of (3 Pic and the impressive accuracy of 
this data minute changes in one's model will result in a 
statistically significant increase in x^- These small in- 
creases in x^ indicate a statistically significant poorer fit 
in the case that one's model is an accurate representation 
of the physical reality of the disk in this case. However, 
in our case it should be noted that even in the best- 
case scenario our model is likely only an approximation 
of the actual physical makeup of the (3 Pic disk; even 

^ For a more de tailed discussion of the r eduction procedure, we 
refer the reader to lGolimowski et al.l I I2006I ) 
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our best-case scenario cannot account for the numerous 
asymmetries in the disk (e.g. Figure [5t . Furthermore, al- 
thoug h the error estimates provided bv lGohmowski et al.l 
(|2006( ) attempt to include both the systematic and sta- 
tistical uncertanty in the ACS data, they could fail to 
account for correlations in the data that could add an 
extra systematic component. For these reasons we feel 
that a simplistic could underestimate the true uncer- 
tainties in our fitted parameters. For these reasons we 
have scaled up the uncertainties in our fitted parameters 
by a factor of 10, as we believe these give a more accurate 
indication of the true uncertainty in our data. 

3. RESULTS 

In the following sections we discuss optimal parame- 
ters for a single and two-disk dust distributions as well as 
several scattering phase functions. Since both the dust 
distribution and the phase function are responsible for 
the scattered light profile, both components are investi- 
gated here and modelled with our simulations. However, 
in order to minimize the number of free parameters, we 
first explore the parameter space using a zodiacal light 
phase function for scattering. In section 3.3 we expand 
on our work by exploring other phase distributions. 

3.1. One Disk vs. Two Disk Model 

To test both the fitting procedure and the goodness- 
of-fit of our two-disk fits, we first constrained the param- 
eters describing the dust distribution via a single disk 
and a zodiacal light scattering profile. The zodiacal light 
profil e has been used in the past (eg. Kalas & Jcwitt 
fl995l) and was chosen in order to minimize the number 
of free parameters. Table [1] lists these results as well as 
those obtained in the past using similar axisymmetric 
models. We do not note any significant disagreements 
and suggest that the minor differences between our and 
past results might be due to several factors including: 1. 
that previou s fitting efforts were p erformed with simpler 
models (eg. iKalas fc JewittjflQQSD . 2. that the previous 
data used for modelling efforts was of poorer quality, 3. 
that various coronographic implementations employed in 
previous observations did not allow the observer to reli- 
ably probe large scale heights, 4. the usage of different 
phase functions and lastly 5. that past modellers applied 
trial and error methods to explore the parameter space 
instead of more sophisticated techniques. 

After we obtained the best single disk (9 parameter) fit, 
we ran our fitting procedure using two disks (18 parame- 
ters). The best- fit results for the one and two-disk models 
are listed in Table [2] along with the associated errors on 
each parameter - the MCMC probability distributions of 
the 18 parameter 606 fit are presented in Figure [3l The 
minimum reduced for the single and two-disk fits are 
1.83 and 1.18 respectively. In addition to the improve- 
ment in the reduced we present other qualitative im- 
provements. In Figure [3] we display the isophotal maps 
of our best-fit models, which shows that the two-disk 
model's isophotes better trace the observed isophotes. 
The isophotal maps also reveal that two disks with ax- 
isymmetric dust distribution and non-isotropic scatter- 
ing profiles can produce asymmetric brightness contours 
if they are inclined to each other and the linc-of-sight 
(LOS) at non-zero angles. Degree of this asymmetry is 
further discussed in section 4.1. 
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Fig. 2. — A plot of the correlation between the radial index (/32) 
and the flare index (72) displaying the 33% (solid thick curve), 68% 
(dotted curve) and the 95% (solid thin curve) credible regions for 
the 18-parameter 606 ACS data. While some correlation betwreen 
these parameters is expected as they both describe light scattered 
at larger radii, the two are still well constrained. 

The total residuals following the removal of the one 
and then two-disk models from the data are presented in 
Figure[5] Aside from producing overall smaller residuals, 
the two-disk model residuals are also symmetric, reveal- 
ing the main warp at (65-85) AU as well as an additional 
structure at (110-150) AU which is similar to the rings 
observed by IWahhai et al.l ()2003f ) . We further discuss 
these structures in section 4.1. In Figure [6] we present 
brightness profile cuts along both the NE and the SW ex- 
tensions of the composite disk, demonstrating once again 
the improvements gained from employing the two-disk 
model. For these reasons we believe the f3 Pic disk is 
much more accurately described by our two axisymmet- 
ric disk model^. 

We also investigate for possible correlations between 
the various parameters of our model. The most signifi- 
cant correlation is between the secondary's outer radial 
index (/32) and the flare index (72). This correlation is 
shown in Figure[2]for the 18-parameter F606W data set. 
Other parameters were not found to be notably corre- 
lated. 

3.2. Chromatic dependence of modeled parameters 

Prior to the iGolimowski et ahl (|2006f ) observations of 
(3 Pic, it was widely accepted that the orbiting dust was 
composed of neutrally scatt ering grains with sizes larger 
than several microns (ex. iLecavelier Pes Etangs et al.l 
119931 ). With the highe r pho tometric precision of the 
ACS. RSolimowski et al.l (|2006l ) analyzed the disk's colors 
with three filters [a narrow B-filter (435 nm), a broadbard 
V (606 nm), and a broadband I (814 nm)] and concluded 
that this assumption was wrong. Complicating matters 
further, their observations showed asymmetries in the 
flux ratios taken at different wavebands about the pro- 
jected major and minor axis of the disk, which imposes 
additional constrains on the models and suggests that a 

^ Bayesian model comparison over the one-disk model while tak- 
ing into accounts the extra degrees of freedom of the two-disk 
model. It was was not applied here due to the clear quantita- 
tive and qualitative improvements of the two-disk model over the 
one-disk model 
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TABLE 1 

Single disk parameter comparison 



Parameter 



Kalas fc Jewitt 



Artymowicz (1998)'^ 



This paper 



Inclination 

Vertical distribution 

Radius of the power law break 

Inner radial index 

Outer radial index 

Flare index 

Scale height 



2° < i < 5° 
0.7 < p < 2.0 
Rm = 100 AU 
N/A 
-3.4 < /3 < -2.8 

1.1 < 7 < 1.6 
0.05 < Ao < 0.10 



P 



a 
7 = 



: 1.3° 

= 0.7 
120 AU 
= 2.0 
: -3.0 
= 0.75 
= 0.055 



i = 2.3° 
p = 0.5 
fl™ = 102 AU 
a = 2.0 
13 = -1.8 
7 = 1.5 
Ao = 0.02 



^ Parameter values converted from IKalas & : Jewittj 19951) to match our model values.'' Personal 
communique. Results obtained by fitting iHeap'etar I (|2000l l HST STIS images. ° Fits obtained by 
trial and error. 



TABLE 2 

9, 18 AND 20 PARAMETER FITS 



Parameter 9 parameter 9 parameter 606 18 parameter 18 parameter 606 20 parameter 20 parameter 606 

606 minimum MCMC best-fit 606 minimum MCMC best-fit 606 minimum MCMC best-fit 





1.0390x10^ 


1.03941°;*^^ X10« 


6.664x10^ 


Reduced x^ 

n n 

pi 

01 (°) 

ai 

pi 

71 


1.8451 
2.261 
0.455 
1.015 

2.04 
-1.84 

1.51 


1.84a»_o oo4g 

9 r,cq + 0.194 

U.400_o 014 
l.UJO_o 
2 05+° " 

, 04-1-0.03 

-l-84_o 02 
1 c:i -1-0.06 
l-''l-0.02 


1.1833 
0.11 
0.84 

0.506 
2.73 

-2.89 
1.05 


i2 n 


n/a 


n/a 


6.312 


P2 


n/a 


n/a 


0.464 


02 (°) 


n/a 


n/a 


3.871 


02 


n/a 


n/a 


1.6 


132 


n/a 


n/a 


-1.16 


72 


n/a 


n/a 


0.84 


Rml (AU) 


102 


102^1 


112.7 


Rm2 (AU) 

VI 


n/a 
3.96x10-1° 


n/a 

4.021°:°^ X 10-10 


71 

5.60x10-10 


m 


n/a 


n/a 


2.0x10-10 


Aoi 


0.018 


U.Ol»_o 002 


0.0489 


Ao2 


n/a 


n/a 


0.051 


9l 


n/a 


n/a 


n/a 


92 


n/a 


n/a 


n/a 



'■664_o.oo2 xlO 



-1-0.0010 
0.0003 
-1-0.06 
0.06 
-t-0.01 
0.02 
-1-0.028 
0.006 
-1-0.12 
-0.08 
,-1-0.07 
'-0.04 
-t-0.04 
0.03 
(. O10-I-0.033 
D.01Z_Q 010 
465+0 024 
U.4DO_o 008 
n 071 -1-0.003 
■J-o( l_o.oo3 
K+0-3 
-0.3 
-1-0.02 
-0.03 
^-t-0.05 
-0.02 

ii2.7i;-o 

1+3 ■ 



1.1834 
0.12 
0.85 
0.505 
2.73 
-2 
1.05 



1.6: 
-1.16 
0.84"' 



7V 



-0.06 
+0.2 
0.1 

-1-0.0010 
0.0012 
-1-0.012 
0.005 



2.0Io:i xlO 
0.0492 
0.052 
n/a 
n/a 



6.6x10^ 
1.17 
0.12 
0.746 
0.48 
4.8 
-2.8 
1.16 
2.6 
1.1 
5.3 
9 

-0.25 
1.59 
112 
63.5 

5.63x10- 
2.41x10- 
0.0429 
0.056 
0.637 
0.849 



6.6j:o2 xlO^ 
-I ,7+0.02 

l-l'-0.04 
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^ As discussed in j|2.4.1l thc formal errors have been scaled up by a factor of 10. 



true description of the dust distribution will be incom- 
plete until wavelength-dependent models are developed. 

Unfortunately, our multi-parametric dust distribution 
model is not wavelength-dependent and neither are the 
phase functions we use'^. In addition, any modifications 
would introduce additional parameters and thus slow 
down the search through the parameter space. Consider- 
ing these points and the already low reduced we have 
been able to obtain from our wavelength-independent 
model, we chose not to increase the complexity of the 
model. 

There are some benefits in our fitting the same model 
to data observed at different wavelengths. Fitting to im- 
ages of the disk at different wavelengths and comparing 

3 zodiacal light function is derived from particles scattering sun- 
light which are much larger than the wavelengths they scatter and 
thus one can think of it as color-neutral scattering function 



the derived parameters can reveal whether the scatter- 
ing in the (3 Pic disk is wavelength dependent and even 
put some constrains on the optical properties of dust. 
Conversely, assuming the scattering in the j3 Pic disk to 
be color-neutral and comparing the results from fits us- 
ing different wavelengths, we can obtain an independent 
estimate of the robustness of our fit. 

Table [3] lists the best fits obtained for the three wave- 
bands. Overall there is a good agreement between the 
parameters, especially for the values describing the pri- 
mary. Several of the parameters describing the secondary 
[the inner radial index a2^ radius of the power law break 
Rm2 and the scale height index ^o] show a much wider 
range than those corresponding to the primary disk. For 
instance the apparent increase of a2 in the shorter wave- 
length bands might suggest that there is substantial red- 
dening in the inner part of the secondary. However, we 
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Fig. 3. — The marginalized likelihood probability distributions of our 18-parameter 606 MCMC chains for each of the fitted parameters. 
The best- fit value (the peak of the distribution) and the 68% credible regions are given by the solid and dashed vertical lines. 



also note that any errors in the subtraction near the 
coronograhic mask would effect a2 the most and thus 
we believe that this result should be viewed with some 
skepticism. 

We also use the models produced for the three filters to 
constrain the inclination of the secondary to the primary. 
The fits for the locations of the spines of both disks (as 
well as the composite) are shown if Figure [71 Table |4] 
lists the corresponding inclinations of the NE and SW 
extensions measured counterclockwise to the major pro- 
jection axis. The value for the angular separation is of 
particular interest since this value has been previously 
reported. Our result, 3.2 ± 1.3°, is cons i stent with the 
~5° value reported bv lGolimowski et"aLl ()2006[ ). 

3.3. Phase Function 

A topic of importance that has largely been ignored 
until now is the appropriateness of different phase func- 
tions to model the scattering of light in the j3 Pictoris 
disk. Since an observed brightness contour is a function 
of both the distribution of scatterers as well as the phase 



function that describes how these particles scatter light, 
a discussion of the latter is necessary. The zodiacal scat- 
tering profile (Fig. [T]) was initially chosen for two rea- 
sons: 1. previou s successful modeling with this function 
(eg. iKalas fc Je witt 1995) and 2. its empirical formula 
does not require any additional parameters. As it was 
already noted, using this profile has allowed us to obtain 
a very low reduced = 1-18. However, the question re- 
mains: how valid is this choice? One way to answer this 
question is to employ other phase functions. In section 
2.2 we introduced the Henyey-Greenstein function along 
with the zodiacal scattering profile. The former is par- 
ticularly interesting because it can be modified through 
the scattering asymmetry parameter g to produce many 
profiles, several of which are shown in Fig. [TJ To obtain 
the best fit, we add the scattering parameter to the 18 
parameters describing dust distribution with two disks 
and run the MCMC algorithm until the different chains 
converge. Intriguingly, the 20 parameter fit does not of- 
fer noticable improvements in the reduced ('^1.2 for 
both 18 and 20 parameters). What's more, the values 
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Fig. 4. — Isophotal map of the hght scattered by the dust orbiting /3 Pic (dotted regions). Also shown are the single (dashed line) and 
the two-disk (solid line) model scattered light profiles. The brightness isophotes are quoted relative to /3Pic, with the isophotes decreasing 
from 1 X 10"'' L^iiak/ Lfjpia with interval spacing equal to 0.5 in log base 10 units. Prom the overall low reduced = 1.84 of a single 
disk fit, one expects a reasonable fit for the observations. This is evident in the isophotal maps. The improvements offered by the two-disk 
model are shown in several regions: along the spine (A), forward scattering areas originating from larger scale heights (B) and the overall 
improvement in tracing the scattered light isophotes (C). 



of the scattering asymmetry profile favored by the pa- 
rameter search {gi = 0.64 and g2 = 0.85 for the two 
disks) produce a forward scattering profile that largely 
resembles the zodiacal scattering function; furthermore 
the uncertainty in gi = 0.64 and g2 = 0.85 are small (see 
Table [U . The disagreement lies in the slightly higher 
amount of back scatter produced by the latter, although 
this is not significant since most of the light comes from 
forward scattering. In Table [5] we present the results of 
our analysis of the two-disk run that allowed for varying 
scattering asymmetry parameters. 

4. DISCUSSION 

4.1. Morphology of the Composite Disk and its 
Components 

The overall shape of the brightness contours produced 
by the two-disk model follows that of the light scattered 
by f3 Pic's dust. This is evident from both the isophotal 
maps (Figured]) and low ~ 1-2 fits for the F606W 
filter images. Similar maps and reduced are recorded 
in other filters as well. The brightness contours of the 
two components (the primary and the secondary) is pre- 
sented in Figure [5] and the relative contribution of each 
component to the total scattered light is shown in Figure 



m 

The most striking feature of these fits is that the sec- 
ondary contributes more light at larger radii and scale 
heights than the primary. The immediate implication is 
that the dust content of the secondary disk is not lo- 
calized to the inner regions of the system and that it 
extends at least o u t to ^ 250 AU. We noted earlier that 
iGolimowski et al.l ()2006( ) found a steeper power law fit 
along the spine of the primary beyond ~ 120 AU than 
the secondary. Their result suggests that the secondary 
will extend further out, at least radially, unless the fit 
is not appropriate beyond ~ 250 AU and is in line with 
our conclusions. It is also interesting that to describe the 
brightness profile along the midplane of the composite, 
IGolimowski et al.l (|20C)6l use 4 power laws with power 
law breaks at 69, 117 and 193 AU. Our two-disk model 
offers an explanation for the locations of the breaks. The 
69 and 117 AU breaks correspond closely to the midplane 
profile breaks of the secondary (71 AU) and the primary 
(113 AU). The last power law break according to our 
model is a result of primary's dust density falling to lev- 
els equal to the secondary's density along the line of sight 
and the contour transitioning to follow the shape of the 
secondary. 
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Fig. 5. — The residual of the 606 ACS data following the removal of the scattered light produced by the single (top panel) and the 2 disk 
(bottom panel) models, presented relative to that of the composite disk. The second residual is not only cleaner, but is also symmetric 
following a 180° rotation. The identified ring-like structures (A and A', B and B') are similar to those identified bv iWahhai et al.l (120031) . 
albeit we are not certain whether their origin is real or failure of the model to properly describe the dust distribution. 



In order to test for possible disagreements with previ- 
ous ground base obs ervations of the disk which extend 
out to 800AU (e.g. iKalas fc Je^dtl Il995( l and do not 
show two separate components, we extend our model to 
the same radius. Due to the large scale height of the 
secondary disk in our model, we find that the composite 
still appears as a single disk at large radial separations 
(out to '^SOOAU). Furthermore, the primary makes a sig- 
nificant contribution to the total light profile along the 
spine, having the effect of aligning the composite's spine 
with its own. Thus, the composite does not show an 
inclination to the primary at 800 AU. 

The impact of such an extended secondary on dust 
dynamics in the system is discussed in section 4.3. 

Figure [3] also reveals that the brightness isophotes of 
the two-disk (composite) model are spaced more widely 
along the northwestern semiminor axis, as can be seen 
from the greater extension of the i s ophot es in the north- 
west direction. iGolimowski et ahl (120061 ) attribute this 
spacing to the northwest side of the disk being tipped 
nearer to earth. All of our two-disk models agree that 
the spacing is a result of the secondary's inclination to 
the line of sight (ia 6.0 ± 1.0° vs h = 0.0 ± 0.1° for 
the primary which is almost entirely edge on). This is 
highlighted in Figure [8] which maps the isophotes of the 
scattered light of the primary and the secondary showing 



that the secondary's isophotes are more widely spaced in 
the NW than the SE, where as the primary fails to shows 
this asymmetry. 

It is also interesting to compare brightness contours 
of the primary and the secondary at different wave- 
lengths. In Figure [10] we present color ratios of the 
F814W and F435W filters f or the primary and the sec- 
ondary. IGolimowski et ahl (|2006D noted that the com- 
posite disk ratio shows asymmetry along both the ma- 
jor and the minor projection axis. The later of these 
was attributed to dust grains scattering blue light more 
isotropically than red light. Our model re-creates this 
asymmetry (Figure llOH op panel) . The ratios of the pri- 
mary at different filters (Figure [TOl-middle panel) do not 
show this asymmetry, which is expected from an axisym- 
metric edge-on disk. Our models suggest that it is the 
secondary which is responsible for the mentioned minor- 
axis asymmetry (which is expected considering its much 
larger inclination to the line-of-sight as well as the pri- 
mary, which was used to define the major and the minor 
axis). 

Lastly, we comment on the residuals obtained from 
the subtraction of the composite from the original im- 
age (Figure HI). As has already been noted, the residual 
is syrn metric resembling rings observed by Wahhai et al] 
(|2003f ) in the mid-IR. The radii of these remnants are 
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Fig. 6. — Scattered light profile along the sciuiiiiajor projection axis (the x-axis of Fig. |4}. Top panels show the northeast and the 
southwest extensions (dashed lined) along with profiles produced by single disk model (thick solid line) and the residual (thin solid line) of 
these two. Bottom panels show the same for the two-disk model. Improvements in the fit can be seen by comparing the two-disk fit with 
the 1 disk fit. The dotted vertical line at 30 AU represents the extent of the coronographic mask. 



(65 - 85) AU and (110 - 150) AU, with the inner ring- 
Hke structure having a similar rad ius and orienta t ion as 
the most outer ring detected by IWahhaj et al.l ()2003l ) 
(R=82 ± 2 AU, i =-2 ± 2°) and coin ciding with the HST 
warp detected bv lHeap et al.l (|2000D . Since this structure 
is detected in the visual and the mid-IR regimes, we do 
not doubt its existence. Furthermore, if the inner "ring" 
is indeed a structure separate from the two disks, it is 
not surprising that our simple axisymmetric descriptions 
fail to account for its presence. Unlike the inner "ring" , 
the outer structure has not been previously observed. 
This raises a possibility that it is a byproduct of the fit- 
ting procedure. For example, the two-disk fit might be 
overcompensating for the inner warp by slightly inclining 
the secondary disk and increasing the overlapping region. 
Our models fail to id entify the 52 AU ring, con firming the 
absence reported bv lGohmowski et"aLl ()2006f ). 



4.2. Asymmetries 

High resol ution images of du s t surr ounding /3 Pictoris 
obtained by iGolimowski et all (|2006D , with their excel- 
lent morphological detail and photometric accuracy, im- 
pose significant constraints on models. In addition, ratios 
of ACS images before and after deconvolution show that 
dust around /3-pic is distributed in two separate disks in- 
clined at ^ 5° to each other. A good model should there- 



fore be able to reproduce the morphology with two nearly 
edge-on disks inclinded by a small angle. As shown al- 
ready, the low reduced and the brightness isophotes 
shown in Figure[4]demonstrate the success of the MCMC 
algorithm we have employed to sweep through the appli- 
cable parameter space. The small inclination of the pri- 
mary to the secondary (3.2 ± 1.3°) and the line of sight 
(ii — 0.0 ± 0.1°) adds additional credibility to our fits. 
However, no model of (3 Pic's dust can be said to be suc- 
cessful unless it recreates at least several asymmetries 
mentioned in the literature, such as the wing tilt and 
the butterfly asymmetry. Below we discuss the known 
asymmetries and the degree to which our model recre- 
ates them. 

Butte rfly Asymmetry: T he asymmetry first iden- 
tified bv lKalas fc JewittI (l995D refers to the asymmetric 
curvature of the isophotes across the spine of the compos- 
ite disk and inversion of this asymmetry across the minor 
axis. Figure m confirms that the two-disk model does a 
good job of re-creating this feature. Also, since the pri- 
mary disk is edge-on and nearly identically aligned with 
the major-projection axis, it is the extended secondary 
that is responsible for this feature^ 

Wing-Tilt Asymmetry: iKalas k Jewit^ (|1995l ) 
noted that the two extensions of the composite disk are 
inclined to each other by 1.3° and that the linear fits for 
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TABLE 3 

Chromatic variations 



Parameter 606 606 MCMC 435 435 MCMC 814 x^ 814 MCMC 

minimum best-fit minimum best-fit minimum best-fit 
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^ As discussed in ^2.4. li the formal errors have been scaled up by a factor of 10, as the uncertainties presented here are believed 
to be a more accurate determination of the uncertainty in the fit parameters. 



TABLE 4 

Relative Position Angles of the Primary, Secondary 
AND THE Composite Disk 



Extension 


F435W 


F606W 


F814W 


Average 


Primary NE 


0.454 


0.478 


0.413 


0.448 


± 


0.033 


Primary SW 


0.460 


0.519 


0.430 


0.470 


± 


0.045 


Secondary NE 


2.567 


2.400 


2.929 


2.632 


± 


0.270 


Secondary SW 


5.342 


5.006 


4.003 


4.784 


± 


0.697 


Composite NE 


0.643 


0.565 


0.446 


0.551 


± 


0.099 


Composite SW 


0.880 


0.782 


0.707 


0.789 


± 


0.086 



All values given in degrees and measured counterclockwise 
from the major projectio n axis. The fits are obtained from 
regions used by Goiimow ski et al.l Hooei), namely 80-250 AU 
for the primary and the composite and 80-130 AU for the sec- 
ondary. 



their midplanes do not intersect at (3 Pic. This asymme- 
try -was named the "wing-tilt asymmetry and is attributed 
to for-ward s cattering from an i ncHne d disk. It is also 
present in the lGoHmowski et all (|2006[ ) images, where the 
inchnation of the primary's components is determined to 
be 0.9° while that of the secondary's is 0.3°. Also, they 
note that the primary's linear fits intersect at a point 
almost coincident with /? Pic, which they note as a point 
of disagreement. 

Our fits also reproduce this asymmetry (Table S]). The 
inclination of the primary's extensions 0.1°), however, 
is mu ch less than that reported by iGolimowski et al"] 
(|2006f ). but is consistent with a disk that is edge-on as 
our model suggests the primary is. The linear fits of 
the disk's spine are nearly centered on the st ar (Figure 
[7]) in agreement with IGolimowski et ahl (|2006) . The sec- 
ondary's extensions on the other hand arc inclined at 
^ 2.2° due to the inclination of this component to the 



line of sight {12 — 6.0 ± 1.0°). Their fits also inter- 
sect close to /? Pictoris. The composite is inclined at 
~ 0.3°, which is less than the inclina t inon of the com- 
posite determined by iKalas fc JewittI (|1995( ). However, 
the two extensions do not meet at the star, with the 
NE component intersecting the SW component ^ 50 AU 
away from (3 Pic along the SW extension. The fact 
that the model primary and composite do not agree on 
the point of intersection is o ffered as a poss i ble re ason 
for the disagreement bet ween iKalas fc JewittI (|1995l ) and 
IGolimowski et al.| ()2006D . Our model suggests that it is 
the combination of the primary and the secondary that 
together creates the ment ioned offset. 

Width Asymmetry: IGolimowski et ah! (^006) find 
that the FWHM and the fuU-width-O.l-max (FWO.IM) 
of the two extensions of the composite diverge beyond 
~ 190 AU and ^ 150 AU respectively. This asymmetry, 
known as the width asymmetry, is present in our fits as 
well, although to a smaller degree. Figure [TT] shows the 
width of the NE and SW components of the composite. 
It is evident that the SW component's FWHM is slightly 
wider than that of the NE component, although not by 
a significant amount. Out to ^ 200 AU the two exten- 
sions F WHM and FWO.IM almos t overlap, as is the case 
with the IGolimowski et al.l (|2006( ) measurements, but af- 
ter this point our mod el disagrees significantly as the 
IGolimowski et aLl ()2006f ) measurements show much more 
asymmetric behaviour. The FWO.IM shows more asym- 
me tric behaviour than the FWHM, as is the case with 
the lGohmowski et all (|2006l ) measurements. 

Surface Brightness Asymmetry: Differences in the 
brightness profile along the spine of the NE and the SW 
extensions of the composite are presented in Figure 1121 
Although our model shows asymmetry between the pro- 
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Fig. 7. — Spines of the composite, primary and secondary disks 
(dashed lines) determined by fitting gaussians to the vertical pro- 
files at each step along the semimajor axis. The solid lines represent 
linear fits to each extension of the disk, which we use to determine 
the inclination values of each component presented in Table [4] The 
solid lines are enlongated into the opposite extension for visibility. 

files of the two extensions, it does this to a much lesser 
degree than the scattered light. 

Radial Extent Asymmetry: Although the radial 
extent asymmetr y could result from th e surface bright- 
ness asymmetry, iKalas fc JewittI ()1995() classify it sepa- 
rately to draw attention to the possibility that the SW 
extension is physically truncated. This asymmetry, how- 
ever, becomes evident only beyond ^ 400 AU and thus 
we are unable to comment on its true nature. 

4.3. Dust Disk Dynamics and Origin 

In previous sections we presented a case for two ax- 
isymmetric disks around (3 Pic. Here we explore the 
dynamical consequences of such a pair. A quantity of 
interest is the mean collision time since it reveals the 
timescales on which the disks are stable, and also on 
which they have to be replenished. It can be estimated 
most easily by using our knowledge about the total stel- 
lar flux absorbed by dust and re-radiated in infrared, 
quantified as the ratio of IR exc ess to star's lumino sity 
fd = LdlL^ = 2.4 X 10-3 fe.g.. IArtvm(wici1fl997l) . If 
the absorbing area of grains in a small volume of space 
is denoted as dAgrain , then we can normalize our model- 



dA 



grain 



kT{r)2TTrdr 

47^7-2 



kT{r)dr 
r 

(10) 

We obtain the normalization constants k for each disk 
separately, utilizing the ratio of scattered light in the 
two disks {Ld,sec/ Ld,prim = 1-23 from our models) and 
the fact that fd is due to the sum of the two disks. The 
normalization factors fcprimary and fcgocondary for the two 
disks can then be used separately or as a sum, to discuss 
the behavior of individual disks or their sum. In fact, 
as there is a large degree of overlap between the disks, 
the second perspective is more rea listic. We estimate the 
mean collision time (tc) following (|Artvmowicz IITqQTI ) as 



tr 



12r(r) 



(11) 



where P{r) is the Keplerian period of a particle at radius 
r. The factor 12 in denominator accounts for the motion 
parallel to the disk and the fact that the cross section 
for collision between two particles of comparable radius 
is up to 4 times bigger than the geometrical cross section 
of either one. 

In Figure [T3] we present the conservative number of 
orbits (assuming the disks completely overlap) a particle 
can make before colliding with another particle as a func- 
tion of distance, and the corresponding mean coUisional 
time scale. 

Since the lifetime of dust in the densest parts of the 
disks is of order 10'' to 10^ yr, and hence much less 
than that of the host star (which is at least 20 Myr 
old), our calculations suggest both disks have planetesi- 
mals embedded in them capable of replenishing the lost 
dust. In other words, none of the disks represent primor- 
dial protoplanetary disks. Rather, they are replenished 
disk(s) owing their existence to a substantial mass reser- 
voir, equivalent to a plan etary system comparable to our 
own (jArtvmowicz |[T997l) . 

Geometrically, the two disks largely overlap and some 
regions contain comparable areas of dust. Physically, 
this implies mutual influence through collisions. Since 
the optical depths of the disks are much smaller than 
unity, particles in the two disks can coexist for as many 
as 20-10000 of their orbits. Similarly, the smaller optical 
depth of the larger parent bodies assures their stability 
for a large number of orbits in the overlapping disks. The 
inferred presence of two disks suggests that (3 Pic has two 
distinct populations of parent bodies (planetesimals) in 
a slightly inclined configuration. 

As mentioned in the introduction, the characterization 
of the distribution of dust around a star is but a first 
step toward the derivation of the mass and spatial distri- 
bution of larger, unseen bodies (planetesimals), as well 
as any perturbing planets; other effects such as dust-dust 
and any gas-dust interaction in the presence of stellar ra- 
diation pressure can spatially displace the fine dust from 
the place of their origin (the parent bodies). However, 
some implications of our models are already apparent. 

The fact that the secondary disk contains a significant 
fraction of the dust in the system means that its ori- 
gin cannot be ascribed to a minor event such as a single 
comet or an asteroid disruption. In order to produce 
comparable amounts of dust in the two disks, they must 



12 



Ahmic et al. 



5 



•a 

o 



■200 




-200 



-50 50 

Projected horizontal distance (AU) 



Fig. 8. — Isophotal maps of light scattered by the primary (top panel) and the secondary (bottom panel). Solid lines are the isophotes 
of the models while the grey regions underneath represent residuals following the removal of the secondary (top panel) and the primary 
(bottom panel) fits. The brightness isophotes axe quoted relative to jBPic, with the isophotes decreasing from 1 x lO"'^ L^igi^/L^pi^ in the 
case of the primary and from 3.162 x 10"'' L^i^k/Lgpi^ in the case of the secondary. In both panels the spax;ing of the isophotes is equal 
to 0.5 in log base 10 units. The empty circular region in the centre of the image represents the coronographic spot covering this area. The 
most important features of these fits include 1 . extended secondary which contributes more light than the primary at larger radii and scale 
heights, 2. more widely spaced isophotes to the NW than SE of the secondary demonstrating signficant inclination and 3. nearly edge-on 
(il = 0.0) primary. 



either have a common massive reservoir of parent bodies 
to draw from or two such reservoirs. Only comprehensive 
modelHng based on high-quahty data can clarify the situ- 
ation since the intermediate-size bodies such as km-sized 
planetcsimals are not amenable to direct detection, pos- 
sessing neither enough gravitational influence nor large 
emitting surface area. 

The vertical profile we obtained is p ~ e"*^!^!/^'""^^ , 
with p ~ pi Ri 0.83 for the first disk and p = p2 ~ 0.46 for 
the second. Thus, neither disk has a vertical scattering 
profile resembling a disk with a vertical isothermal struc- 
ture = 2, or Gaussian vertical profile). The two disks, 
especially the second one, have super-exponential drop- 
offs (p < 1) and therefore extended, slowly diminish- 
ing outer wings of the vertical structure, with a sharply 
pointed inner core. This could be a sign that the inner 
core of the profile, referred to as the spine in this paper, 
is composed of bigger particles than the profile wings; 
the mass separation might arise due to energy equiparti- 
tion leading to a smaller velocity dispersion of the larger 
particles in the disk. In addition, a more fanned out sec- 
ondary debris disk could be a consequence of more en- 



ergetic collisions of its particles, perhaps indicating their 
generally smaller sizes and larger radiative modification 
of orbits. 

An expectation based on the the previous dynamical 
studies of the mobility of dust is that the bulk of a disk 
or a ring of planetcsimals will be located interior to the 
dust disk. This follows from the outward action of ra- 
diation pressure exerted on dust. If so, we expect that 
most planetcsimals will be found at a distances less than 
50-100 AU from (3 Pic. At such distances, planetcsimals 
can be perturbed by the proposed 8- Jupiter mass planet 
at a mean distance of 8 AU and/or other unknown plan- 
etary masses. Perturbing planets, in any case, must be 
located inside the cross-over radius {Rm2 < 70 AU) of 
the secondary disk. 

We stress here the possibility, indeed necessity, of mul- 
tiple planets in [3 Pic. A planetary system born with 
only one major planet would tend to stay in alignment 
with the orbital plane of this planet. In such a case, 
we would not expect to see two separate inclined disks. 
What could be the origin of the inclined orbit of an ad- 
ditional planet? Although specific inclination pumping 
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Fig. 9. — Ratio of the projected light scattered by the two disks, 
demonstrating that along the midplane light from the primary 
dominates, while although at larger scale heights the secondary 
becomes dominant. The solid line outlines the region inside of 
which the primary's light dominates. The dashed line inside the 
sold lines denotes the isophote where the primary's contribution is 
twice that of the secondary. The dashed lines outside the region 
represent isophotes where the secondary's contribution is twice and 
ten times greater than the primary's, the latter being the isophote 
at larger scale height. 

mechanisms have yet to be quantified, this may provide 
a natural explanation to the inclination of planetesimals 
and dust from the second disk. This could result in the 
planet (s) themselves being on slightly inclined, mutually 
interacting, orbits thus contributing to the survival of 
the two inclined disks. 

Future observations of the inner regions of the disks 
would be instrumental in establishing the presence of 
the putative multiple planets via their direct influence on 
dust. At some point, the assumption made in the present 
work that the disks are intrinsically axisymmetric, if mis- 
aligned, must prove inadequate. With sufficient observa- 
tional data, our assumption of the same optical proper- 
ties of dust in two disks may also be relaxed. Further- 
more, a future comprehensive model of the system must 
include not only coUisional dust dynamics and radiation, 
but also the coUisional cascades spanning a large size 
range from planetesimals to micron-sized dust, and pos- 
sibly more unusual ingredients, such as dust avalanches 
involving sub- micron dust. 

5. CONCLUSIONS 
Below we summarize the results of this work. 

1. We modeled the dust distribution around (3 Pic 
with two axisymmetric disks. A MCMC op- 
timization method was used to reject disk models 
and phase functions that produce brightness con- 
tours differe nt from t hose obtained with the ACS 
bv lGolimowski et al.l ([2006). We find a small quan- 
titative improvement in the reduced (from 1.8 to 
1.2) of the two-disk fit compared to the single disk 
fit; the two-disk fit, however, displays significant 
improvements in its ability to trace the isophotes 
of the composite. The parameters listed in Table 
[3] along with a zodiacal light scattering profile pro- 
duce the best fits. 
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Fig. 10. — Ratios of light scattered by the composite (top panel), 
primary (middle panel) and secondary (bottom panel) disks at 435 
and 814 nm. The red region shows the highest relative contribution 
in the 814 nm relative to the 435 nm filter. The ratios are presented 
in log base 2 units. From the presented ratios, we conclude that 
the color asymmetries along the projected major and minor axis 
are due to an inclined secondary disk. 

2. Dust distribution in the primary disk is largely 
in agreement with conclu sions from previou s mod- 
eling attempts (eg. lArtvmowicz et al] 119891 : 
iKalas fc Je\^ll995f) . The modelled primary disk 
has a near-zero inclination to the line of sight. 
Therefore, the up-down asymmetry of the whole 
disk is not caused by the primary. The surface 
density of dust increases as ai = 2. 73-th power 
of radius inside a cross-over radius Rmi =113 AU 
and decreases at large radii with power-law index 
Pi = -2.89 (values from F606W filter fit). This 
behavior is similar to the one known from previous 
modelling, although the inner power law is steeper. 
The vertical disk profile has a nearly exponentiall 
fall-off with vertical distance {pi = 0.85). The 
scattered light contours from the primary disk are 
shown in Figure [H 

3. The secondary disk model presents the first quan- 
titative measurement of the dust distribution in 
this component. It suggests that the secondary 
disk, with a scale-height much larger than the pri- 
mary, is not radially limited to the inner system. 
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Fig. 11.— FWHM (lower pair) and the FWO.IM (upper pair) 
of our model for the NE (solid lines) and SW (dashed lines). The 
wider SW extension demonstra tes that our model rep roduces the 
width asymmetry first noted bv lKalas fc JewittI (119951) . 
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Fig. 12. — The solid line represents the ratio of the brightness 
profile along the midplane of the SW extension divided by the 
corresponding brightness profile of the NE extension of our best-fit 
18-parameter model to the 606 ACS data. The same ratio is shown 
for the 606 ACS scattered light data (open circles). All open circles 
represent a 3 pixel average taken along the minor projection axis 
(the y-axis of Figure 4). The model profile along the spines of 
the two extensions shows asymmetric behaviour, but the degree of 
this asymmetry is far less than that observed from the scattered 
light. The dotted vertical line at 30 AU represents the extent of 
the coronographic mask. 

but that it extends at least to 250 AU, and pos- 
sibly further out. The secondary disk model has 
a larger scattering area than the primary at large 
distances and scale heights, which follows from the 
less-steep power laws of the radial surface density 
{(32 = —1.16 > —2.89 at distances much larger than 
the cross-over Rjn2 = 71 AU) and the vertical den- 
sity {p2 = 0.46 < 0.85). The ratio of the total 
scattered light (secondary vs. primary) is roughly 
equal to 1.23. Since we do not model the dust size 
distribution, we cannot compute the dust mass ra- 
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Fig. 13. — The optical depths and mean coUisional times for 
the two disks (assuming they are completely physically separated; 
dashed lines for the primary, dotted lins for the secondary) as well 
as the conservative estimate (assuming that the two disks overlap 
completely; solid lines). Due to the uncertainty in our model of the 
inner radial index (02), the coUisional times for the inner regions 
of the secondary disk (<~ 60 AU) are less reliable. 

tio of the two components. It is however plausible 
that the secondary disk is as substantial as the pri- 
mary. The larger extent and especially the larger 
scattering area of the secondary disk are surpris- 
ing and raise many questions about the origin and 
dynamics of such a pair of disks. The scattered 
light contours from the secondary are presented in 
Figure [H 

4. Several of t he asymmetries reported by 

iKalas fc JewittI (|1995[ ) can be reproduced with 
our model. The butterfly asymmetry is explained 
by the presence of an extended secondary disk 
inclined by 3.2 ± 1.3° to the primary. We find that 
it is the secondary that is inclined to the line of 
sight and therefore responsible for the wing-tilt 
asymmetry observed in the composite image. The 
width as ymmetry is present , but n ot to the degree 
seen in iGolimowski et al.l (|2006D data beyond 
~ 200 AU. The brightness asymmetry is also 
noticable both across the major and the minor 
projection axis. The difference across both the 
axes is attributed to a secondary disk inclined 
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by 12 = 6.0 ± 1.0° to the line of sight and to the 
primary which is nearly aligned with the major 
projection axis. 

5. In order to explore dust dynamics in a pair of disks 
described by the parameters obtained here, we an- 
alytically constrain the dust replenishment times 
to be 10^ yr at ^100 AU. Such small timcscales 
show that both disks have to be continuously re- 
plenished and contain planetesimals that, as even- 
tual mass reservoir, are capable of producing the 
second-generation dust. The nature and origin 
of the inclined disks and planetesimals in (3 Pic 
should be explored via dynamical models of spa- 
tially resolved coUisional cascades subject to radi- 
ation pressure effects. 

6. The most plausible conjectural outcome of this 
modeling is that the non-coplanarity of two dust 



Artymowicz, P., Burrows, C, Parcsce, F., 1989, ApJ, 337, 494 
Artymowicz, P., 1997, Aimu. Rev. Earth Planet. Sci., 25, 175 
Aumann, H.H., ct al. 1984, ApJ, 278, L23 
Amnann, H.H., 1985, PASP, 97, 885 

Crifo, F., Vidal-Madjar, A., Lallcmcnt, R., Fcrlet, R., Gcrbaldi, M., 

1997, Astron. Astrophys., 320, 29 
Crifo, F., Vidal-Madjar, A., Lallcment, R., Ferlet, R., Gerbaldi, M., 

1997, Astron. Astrophys., 320, 29 
CroU, B. 2006, PASP, 118, 1351 
Ford, E.B. 2005, AJ, 129, 1706 

Gelman, A., Rubin, D.B. 1992, Statistical Science, 7(4), 457 
Gillct. F., 1986, in Light on Dark Matter, edited by P.P. Israel 

(Reidol Dordrecht), p. 61 
Golimowski, D. A., Durrance, S. T., & Clampin, M. 1993, ApJ, 

411, L41 

Golimowski, D. A., ct al. 2006, AJ, 131, 3109 

Grigoricva, A., Artymowicz, P, & Thcbault, Ph. 2007, A&A, 461, 
537 

Heap, S. R., Lindler, D. J., Lanz, T. M., Cornctt, R. H., Hubeny, 

I., Maran, S. P., & Woodgate, B. 2000, ApJ, 539, 435 
Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70 



disks is due to the secondary dust disk being sup- 
plied by its own disk of planetesimals, which has 
been inclined by the gravity of a planet residing in- 
side the cross-over radii (r < 70 AU). This, in turn, 
would be most naturally explained if there were 
more than one massive planets in the system, in 
orbits inclined by several degrees. Future modeling 
of unseen planetesimals of /3 Pictoris, constrained 
by the new knowledge of dust distribution, may 
clarify this conjecture. 
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